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Abstract 

The debonding of skin/stringer constructions is analyzed using a step-by-step 
simulation of material degradation based on strain softening decohesion elements and 
a ply degradation procedure. Decohesion elements with mixed-mode capability are 
placed at the interface between the skin and the flange to simulate the initiation and 
propagation of the delamination. In addition, the initiation and accumulation of fiber 
failure and matrix damage is modeled using Hashin-type failure criteria and their 
corresponding material degradation schedules. The debonding predictions using 
simplified three-dimensional models correlate well with test results. 


INTRODUCTION 

The fracture process of high performance composite laminates is quite complex, 
involving both intralaminar damage mechanisms (e.g. matrix cracking, fiber fracture) 
and interlaminar damage (delamination). The basic failure mechanisms that can 
induce failure in a laminated composite are: matrix tensile/compressive cracking, fi- 
ber-matrix shearing, fiber tensile failure, fiber microbuckling, and delamination. To 
simulate damage growth accurately, a failure analysis must be able to predict all 
failure modes and apply a corresponding reduction in material stiffness. 

The objective of the proposed paper is to demonstrate a methodology that can cal- 
culate the strength of composite structures by predicting the growth of all major dam- 
age modes, including delamination and ply damage. The effectiveness of the method- 
ology for modeling the debonding of skin/stringer is assessed. The failure of these 
specimens is clearly dominated by delamination. However, other damage mechanisms 
such as transverse matrix cracking are also present and the simulation of their interac- 
tion with the initiation and growth of delamination needs to be included in the analy- 
ses. Delamination is simulated with decohesion elements and ply damage is simulated 
by degrading the elastic properties after a failure criterion has been exceeded. 
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DECOHESION ELEMENTS 


Decohesion elements are well suited to progressive failure analyses where both 
ply damage and delaminations are present. The approach consists of placing interfa- 
cial decohesion elements between composite layers. A decohesion failure criterion 
that combines aspects of strength-based analysis and fracture mechanics is used to 
simulate debonding by softening the element. The constitutive equations for the inter- 
face are phenomenological mechanical relations between the tractions and interfacial 
separations. With increasing interfacial separation, the tractions across the interface 
reach a maximum, decrease, and vanish when complete decohesion occurs. The work 
of normal and tangential separation can be related to the critical values of energy re- 
lease rates. 

The decohesion elements are used to model the interface between sublaminates or 
between two bonded components. An 8-node decohesion element shown in Fig. 1 was 
developed and implemented in the ABAQUS finite element code. The element is a 
zero-thickness volumetric element in which the interpolating shape functions for the 
top and bottom faces are compatible with the kinematics of the elements that are being 
connected to it. Details of the development of this element are presented in Refs. 1 
and 2. 



Constitutive Equations of Decohesion 

A constitutive equation is used to relate the traction <7 to the relative displacement 
8 at the interface. Some strain softening models that have been proposed are shown in 
Fig. 2 and include: linear elastic-perfectly plastic; linear elastic-linear softening; linear 
elastic-progressive softening; linear elastic-regressive softening; and Needleman [3]. 
One characteristic of all softening models is that the cohesive zone can still transfer 
load after the onset of damage (5° in Fig. 2). For pure Mode I, II or III loading, after 
the interfacial normal or shear tractions attain their respective interlaminar tensile or 
shear strengths, the stiffnesses are gradually reduced to zero. The area under the 
traction-relative displacement curves is the corresponding (Mode I, II or III) fracture 
energy. Using the definition of the / integral proposed by Rice [4], and evaluating the 
integral Jon the contour defined by the cohesive zone, the relationship 
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is obtained, where G c is the critical energy release rate for a particular mode and 8' is 
the corresponding relative displacement at failure (8 PP , 8 pro , <5/„, <5v e , or 8 re in Fig. 2). 



- Perfectly plastic (pp) 

- Linear softening (lin) 

- Progressive softening (pro) 

- Regressive softening (reg) 
Needleman (Ne) 


Figure 2. Strain softening constitutive models. 


Bilinear Softening Model 

The linear elastic-linear softening (bilinear) model is the simplest to implement, 
and is most commonly used [1; 5; 6]. The double cantilever beam test shown in Fig. 3 
illustrates the material response. The stress <x represents a normal or tangential trac- 
tion between the top and bottom faces of the element, which is a function of the rela- 
tive normal or tangential displacement 8. Point 1 in the plot and the DCB shown in 
Fig. 3 corresponds to a material location subjected to a low tensile load that is within 
the elastic range. A high initial stiffness K p (penalty) holds the top and bottom faces of 
the interface element together. Point 2 represents the onset of damage. In single-mode 
delamination, the traction at point 2 is equal to the corresponding interlaminar strength 
of the material, o c . As the relative displacement increases, the interface accumulates 
damage and the traction is lower than the strength (point 3). The energy dissipated 
damaging the surface at point 3 is the area of the triangle 0-2-3. The energy released is 
the most important component of the failure criterion described in the next section. If 
the load were to reverse, point 3 would unload to the origin, as shown in the figure. 

The critical value of the energy release rate is attained at point 4. For any relative 
displacement larger than point 4, the interface does not carry any tensile or shear loads 
(point 5). In other words, at point 4 all the available interfacial fracture energy has 
been completely consumed. Note that when modeling delamination with a softening 
response, the delamination tip is not defined explicitly. While the onset of damage 
occurs at point 2 in Fig. 4, the delamination tip could be defined as the point where the 
tractions at the interface are zero, which is point 4. 

The softening response illustrated in Fig. 3 is representative of the tension or the 
shear response but not compression. It is assumed that compression loads do not cause 
delamination or softening, and the effect of compression on damage of the interface 
was neglected in the present work. 


p 



Figure 3. Bilinear constitutive model. 


The problem of contact of the crack faces after failure is addressed by re-applying 
the normal penalty stiffness. The process of reapplying the normal stiffness when 
interpenetration is detected is typical of solution procedures of contact problems using 
penalty methods in a constrained variational formulation. 

Mixed-Mode Delamination Criterion 

Under pure Mode I, 11 or 111 loading, the onset of damage at the interface can be 
determined by comparing the stress components with their respective allowables. 
Under mixed-mode loading, however, damage onset may occur before any of the 
stress components involved reach their respective allowables. The mixed-mode 
criterion proposed by the authors in Ref. 2 assumes delamination initiation can be 
predicted using a failure criterion based on a quadratic interaction between the 
interlaminar stresses. In addition, it is also assumed that the delamination growth is 
governed by a quadratic interaction among the energy release rates corresponding to 
the opening (Mode I) and tangential (Modes II and III) displacements. 

The element stiffness matrix is calculated based on the standard isoparametric 
linear Lagrangian interpolation functions for three-dimensional (8-node) elements. 
The element was implemented into an ABAQUS user- written subroutine (UEL)[7]. 
The details of the formulation can be found in Refs. 1 and 2. 


PLY DAMAGE 

To simulate damage growth accurately, the failure analysis must be able to predict 
the failure modes in each ply and apply the corresponding reduction in material 
stiffnesses. The failure criteria included in the present analyses are those proposed by 
Hashin [8] and are summarized below. 

Matrix Failure in tension and compression occurs due to a combination of 
transverse stresses, O22, O53, and shear stresses, T12, in and t 23. The failure index can 
be defined in terms of these stresses and the strength parameters Y and shear 
allowables S c . Failure occurs when the index exceeds unity. Assuming linear elastic 
response, the failure index has the form: 
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where is the strength perpendicular to the fiber direction in tension, Y c is the 
strength perpendicular to the fiber direction in compression, and S c n, S c i 3 , and S C 23 are 
the in-plane shear and transverse shear strengths, respectively. 

Fiber-matrix failure occurs due to a combination of axial and shear stresses. The 
failure criterion has the form: 
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where X, is the strength along the fiber direction in tension, and X c is the strength 
along the fiber direction in compression. 

Fiber failure in tension or compression occurs independently of other stress 
components. In compression the fiber fails by buckling. The failure criterion has the 
form: 


e f = 
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and e f =^~ for < 7 n >0 (7) 
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To simulate the above failure modes, the elastic properties are made to be linearly 
dependent on three field variables, FV1 through FV3. The first field variable 
represents matrix damage, the second the fiber-matrix shearing, and the third the fiber 
buckling failure. The values of the field variables are set to zero in the undamaged 
state. After the failure index has exceeded 1.0, the associated user-defined field 
variable is set to 1.0. The associated field variable then continues to have the value of 
1.0, even though the stresses may reduce to values lower than the failure of the 
material. This procedure ensures that the damaged material does not “heal”. The 
mechanical properties in the damaged area are reduced appropriately, according to a 
pre-defined degradation schedule. For example, when the matrix failure criterion takes 



the value of 1.0, the transverse shear modulus E y and Poisson ratio Vn are set equal to 
zero. The field variables can be made to transition from 0 (undamaged to 1 (fully 
damaged) instantaneously. 

The finite element implementation of this failure analysis was developed for the 
ABAQUS structural analysis program using the USDFLD user- written subroutine [7], 
ABAQUS calls the USDFLD subroutine at all material integration points of elements 
that have material properties defined in terms of the field variables. The subroutine 
provides access points to a number of variables, such as stresses, strains, material 
orientation, current load step, and material name, all of which can be used to compute 
the field variables. Additional details on the procedure can be found in Ref. 9. 

The procedure for the simulation of ply damage can be applied in conjunction 
with decohesion elements to model delamination growth. The following sections 
describe analyses of skin/stringer debonding in which ply damage and delamination 
growth are modeled simultaneously. 


ANALYSES OF SKIN-STRINGER SPECIMEN 

Extensive testing of the skin-stringer flange debonding problem has been 
performed by Cvitkovich and O’Brien [10; 11], and analyses have been conducted by 
Krueger and O’Brien [12], The configuration of the specimens studied in Refs. 10 and 
12 is shown in Fig. 4. The specimens are 203 mm.-long and 25.4 mm.-wide. Both 
skin and flange are made from IM6/3501-6 graphite/epoxy prepreg tape with a 
nominal ply thickness of 0.188 mm. The skin lay-up consisting of 14 plies is 
[0/45/90/-45/45/-45/0] s and the flange lay-up consisting of 10 plies is [45/90/- 
45/0/90] s . Two loading configurations were considered: tension and three-point 
bending. 
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Figure 4. Skin-stringer specimen and test configurations. 



Reference 12 describes the complexities of this problem, such as the presence of 
several interacting failure modes, a non self-similar crack growth, and crack fronts 
that jump between plies and do not remain at a single interface. References 10 and 12 
also show that delamination seldom initiates or grows between the flange and the skin, 
but rather on a surface one or two plies away from the skin-flange interface. A photo 
micrograph taken from Ref. 10 and shown in Fig. 5 shows these features. 

The moduli and strength properties of the unidirectional graphite/epoxy and the 
the interface reported in Refs. 9 and 12 are shown in Tables 1, 2, and 3, respectively. 



Figure 5. Detail of debond area in a failed specimen [1 1], 


Table 1. Elastic properties of unidirectional IM6/3501-6 graphite/epoxy [12]. 


E n 

E 22 -E 33 

Vl 2 — Vl3 

V23 

G 12 —G 13 

G 2 3 

144.7 GPa 

9.65 GPa 

0.3 

0.45 

5.2 GPa 

3.4 GPa 

Table 2. Strength properties of unidirectional IM6/3501-6 graphite/epoxy [9]. 

X t 

Xc 

Y,[ 12] 

Y c 

S 12~S 13 

S 23 

1450 MPa 

1130 MPa 

61 MPa 

250 MPa 

142 MPa 

68 MPa (est.) 


Table 3. Properties of the interface [12]. 

Gi c G IIc Z t 

$lt 

0.075 N/mm 0.60 N/mm 6 1 MPa 

68 MPa (est.) 


Simplied Model of Skin-Stringer Specimen 

To perform fracture mechanics calculations of skin/stringer debonding, Krueger 
[12] developed detailed two-dimensional models with up to four elements per ply 
thickness in which the discrete matrix cracks and the specific delamination paths ob- 
served in the experiments were modeled. Such level of detail may not be practical for 
design analyses, especially for components in which the crack paths have not been 
determined experimentally beforehand. 

The objective of the present work is to demonstrate a methodology with predictive 
capabilities applicable to larger structural components. Consequently, simplifications 
were performed to ease the development of the model and to reduce the number of de- 


grees of freedom while retaining the ability to predict accurately the debond load and 
the subsequent response. A relatively coarse three-dimensional model with 15,212 
degrees-of-freedom was developed in which the thickness of the skin and the flange 
are represented by only two brick elements each. The delamination was assumed to 
occur at the interface between the skin and the flange, although the experiments 
indicate that it always develops one or two plies above or below the interface. To 
prevent delamination from occurring at both ends of the flange simultaneously, 
symmetry in the model was perturbed by modeling the taper of the flange with a 
refined mesh on one end and a coarser mesh on the other. 

Unsymmetric layered properties were assigned to the three-dimensional elements 
by using the composite solid section option in ABAQUS. This option calculates the 
element’s stiffness by performing a Simpson’s integration through the thickness using 
material points in each of the plies within the element. Approximations were 
necessary in the tapered sections of the flange because planes of integration cannot be 
defined parallel to the plies in the tapered elements, as illustrated in Fig. 6. 


Geometry of Simplified 



^ Decohesion elements 

Figure 6. Detail of flange tip geometry and corresponding FE model. 


Interacting Matrix Cracking and Delamination 

The methodology proposed herein models the simultaneous accumulation of ply 
damage and delamination growth. Furthermore, the load incrementation procedure 
represents the interactions among the different failure modes by taking into account 
the load redistribution that results from any type of damage. However, a continuum 
damage model does not adequately represent the interaction between transverse ma- 
trix cracks and delamination. For instance, when the matrix crack in Fig. 5 reaches the 
45° ply beneath it, the crack turns and initiates a delamination. The effect of matrix 
cracks on the model can be the opposite: a flange element with matrix cracks softens 
and it no longer carries the loads necessary to cause a delamination. To address the 
issue, it was assumed that the decohesion elements outboard of a flange element with 
transverse matrix cracks should be removed. Unfortunately, the current version of 
ABAQUS Standard (V6.2) does not support automatic element erosion. Therefore, it 
was necessary to manually disable the corresponding elements and to restart the 
analysis from the load increment in which the matrix cracks were detected. 


Simulation of Tension Test 


The boundary conditions for the tension test are shown in Fig. 4. The mechanism 
for debonding under axial load is primarily by shear lag. Deformed plots of the finite 
element model immediately before and after flange separation are shown in Fig. 7 
with a displacement magnified 8 times. The figure shows that the refined end of the 
flange debonds and the coarse end does not. A detailed analysis of the results indi- 
cated that the coarse end of the flange also softens, but that the separation of the 
flange at the refined end relieves the stresses at the other end. Interestingly, static tests 
exhibit debond of only one end, while fatigue tests induce debonding at both ends. 

Figure 7 also shows that the debond growth is not symmetric: the debond initiates 
on the left comer of the flange shown on the top of Fig. 7 due to the unsymmetry in- 
troduced by the terminated plies at the flange tapered ends. Cvitkovich [10] also 
commented on the difference in the modes of failure of the two comers. 




Figure 7. Deformed plots of skin/stringer model at applied loads of 15 and 28 kN (8X magnified). 




integration points 


Figure 8. Matrix failures at the tip of the flange. 

The damage initiation mode is by matrix failure at the 45-degree ply of the first 
row of elements, as shown in Fig. 8. It can be observed that the damage is non- 
symmetric due to the 45-degree ply which is longer at the taper, as shown in Fig. 6. 

The axial strains from the two gages shown in Fig. 4 are are compared to the 
present analytical results in Fig. 9. The strains and the failure load predicted by the 
model correlate well with the test results. However, the test results after the debond 
load are misleading because the load was immediately reversed after failure by the test 
operator, even if the specimen could carry a higher load. The predicted load at failure 
is 23.9 kN, which is only 1.3% higher than the highest debond load measured. It is 
interesting to note that the debond load predicted by ignoring all damage modes ex- 
cept delamination is 24.1 kN, which con fir ms that the failure of the skin/stringer 
specimens is dominated by delamination. 

The extensometer measurements for tests C4, C6, C8 and CIO conducted by 
Cvitkovich [10] and the predicted values are compared in Fig. 10. The initiation of 
delamination is marked by the sharp breaks in the extensometer readings. After fail- 
ure, the extensometer measurements are unreliable, as can be observed by comparing 
the different tests. The predicted response correlates well with the experimental re- 
sults. The predicted debond load is only 4% higher than the range for all four tests. 



Figure 9. Strains at two gage locations for test and analysis (tension test). 



Figure 10. Predicted and experimental extensometer measurements for tension test. 


Simulation of Bending Test 

The boundary conditions for the three-point-bending test are shown in Fig. 4. The 
mechanism for the initiation of debonding under axial load is a combination of peel 
and interlaminar shear. The initial damage event consists of matrix cracks in the 90- 
degree ply of the flange taper at a load of 405 N (damage location #1 in Fig. 11). It is 
assumed that the row of elements at the flange tip do not carry loads, so that the first 
row of decohesion elements (row A) is manually eroded, as explained previously. The 


second damage event occurs at 460 N, and it consists of matrix cracks in the 45- 
degree ply (damage location #2 in Fig. 11). At this load, the next two rows of 
decohesion elements are eroded (rows B and C). 

Undeformed and deformed plots of the finite element model are shown in Fig. 12, 
where the displacements have been scaled by 200% for clarity. At the final mid-span 
displacement of 6.5 mm., nearly half of the flange has debonded. 



Figure 11. Sequence of initial damage events at flange tip (load=460 N.) 


Deformed d=6.5 mm. 



Figure 12. Undeformed and deformed plots of skin/stringer model in three-point bending. 


The load-deflection curves for the tests and the simulation are shown in Fig. 13. It 
can be observed that the analysis overpredicts the debond load by 20%. The analysis 
predicts an initial unstable delamination (load reduction) of approximately % of the 
span of the flange, followed by additional marginally stable growth at a nearly 
constant load of 600 N until half of the flange has debonded. In contrast, the test 
specimens exhibited an unstable growth of approximately half of the span of the 
flange. The difference between the experimental observation and the prediction can be 


attributed to the dynamic effects of crack propagation, which are not considered in the 
simulation. 



Figure 13. Predicted and experimental load-deflection curves for bending test. 


CONCLUDING REMARKS 

A methodology that can simulate the accumulation of ply damage and 
delamination growth was developed. Decohesion elements were used to model the 
initiation and growth of delamination, while a more conventional property 
degradation procedure was used to simulate ply damage modes such as matrix 
damage and fiber failures. The methodology was applied to the prediction of the 
debond loads of skin and streinger, a problem clearly dominated by delamination but 
where other damage mechanisms are also active. Two specimen configurations were 
analyzed: a tension test and a three-point-bending test. The predicted debond loads for 
both configurations agree well with the test results. However, it was found that a 
continuum damage model for the transverse matrix cracks does not provide the 
adequate conditions to initiate the delaminations observed in the tests. The difficulty 
was successfully resolved by manually initiating a delamination under the elements 
exhibiting transverse matrix cracks. 
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